#include "ManagedLasoProvider.h"

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::slaeig_
   (integer *n, integer *nband, integer *nl, integer*nr, real *a, real *eigval,
   integer *lde, real *eigvec, real *bound, real *atemp, real *d__, real *vtemp,
   real *eps, real *tmin, real *tmax)
{
  /* System generated locals */
  integer a_dim1, a_offset, eigvec_dim1, eigvec_offset, i__1;
  real r__1, r__2;

  /* Builtin functions */
  //double sqrt(doublereal);

  /* Local variables */
  static real atol;
  static integer nval, i__, m;
  static real artol;
  //extern /* Subroutine */ int slabcm_(integer *, integer *, integer *, 
  // integer *, real *, real *, integer *, real *, real *, real *, 
  // real *, real *, real *, real *);


  /*  THIS IS A SPECIALIZED VERSION OF THE SUBROUTINE BNDEIG TAILORED */
  /*  SPECIFICALLY FOR USE BY THE LASO PACKAGE. */


  /*  LOCAL VARIABLES */


  /*  FUNCTIONS CALLED */


  /*  SUBROUTINES CALLED */

  /*     SLABCM, SLABFC, SLAGER, SCOPY */

  /*  SET PARAMETERS */

  /* Parameter adjustments */
  a_dim1 =  *nband;
  a_offset = 1+a_dim1 * 1;
  a -= a_offset;
  --eigval;
  eigvec_dim1 =  *lde;
  eigvec_offset = 1+eigvec_dim1 * 1;
  eigvec -= eigvec_offset;
  bound -= 3;
  --atemp;
  --d__;
  --vtemp;

  /* Function Body */
  /* Computing MAX */
  r__1 =  *tmax, r__2 =  - (*tmin);
  atol = (real)(*n) **eps * max(r__1, r__2);
  artol = atol / sqrt(*eps);
  nval =  *nr -  *nl + 1;

  /*   CHECK FOR SPECIAL CASE OF N = 1 */

  if (*n != 1)
  {
    goto L30;
  }
  eigval[1] = a[a_dim1 + 1];
  eigvec[eigvec_dim1 + 1] = (float)1.;
  return 0;

  /*   SET UP INITIAL EIGENVALUE BOUNDS */

  L30: m = nval + 1;
  i__1 = m;
  for (i__ = 2; i__ <= i__1; ++i__)
  {
    bound[(i__ << 1) + 1] =  *tmin;
    bound[(i__ << 1) + 2] =  *tmax;
    /* L50: */
  }
  bound[4] =  *tmax;
  bound[(nval + 2 << 1) + 1] =  *tmin;
  if (*nl == 1)
  {
    bound[4] =  *tmin;
  }
  if (*nr ==  *n)
  {
    bound[(nval + 2 << 1) + 1] =  *tmax;
  }

  /* L60: */
  slabcm_(n, nband, nl, nr, &a[a_offset], &eigval[1], lde,
    &eigvec[eigvec_offset], &atol, &artol, &bound[3], &atemp[1], &d__[1],
    &vtemp[1]);
  return 0;
} /* slaeig_ */
